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A simple formula is obtained for coupling electrons in a complex system to the electromagnetic 
field, ft includes the effect of intra-atomic excitations and nuclear motion, and can be applied in, e.g., 
first-principles-based simulations of the coupled dynamics of electrons and nuclei in materials and 
molecules responding to ultrashort laser pulses. Some additional aspects of nonadiabatic dynamical 
simulations are also discussed, including the potential of "reduced Ehrenfest" simulations for treating 
problems where standard Ehrenfest simulations will fail. 
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It is now possible to perform first-principles simu- 
lations of the coupled dynamics of electrons and nu- 
clei with all nuclear coordinates included [H, [H, H 0], 
rather than a subset of nominal reaction coordinates. 
For very large systems, or when many trajectories are 
necessary, it is convenient to use a first-principles- based 
scheme [S H, 0, HI, with a valence-electron Hamilto- 
nian and ion-ion repulsive potential derived from calcu- 
lations using density-functional or other first-principles 
techniques. Here we are mainly concerned with the is- 
sue of how one can efficiently and accurately couple elec- 
trons to the electromagnetic field in such an approach, 
where matrix elements of various operators between lo- 
calized basis functions (or "atomic orbitals" ) can be cal- 
culated from first principles, and then used in large-scale 
calculations for complex systems, such as materials and 
molecules, responding to applied fields, such as ultrashort 
laser pulses [g EE [Til El El HI EE 111, El El El ■ 

Our starting point is, of course, the time-dependent 
Schrodinger equation 

ih^-i>(x,t) = Hi>(x,t) (1) 
at 

H = —(-iKV--A(x,tj) +U , q=-e.(2) 
2m V c / 

Some time ago, Graf and Vo gl l2Cl|| obtained a result, 
used in Refs. |ia.[li.[T5l.[l^[T^1l8lll9j|. which is the time- 
dependent version of the Peierls substitution: If H is the 
Hamiltonian matrix in a localized basis with no applied 
field, 



H (£', I) = / d 3 x (x - X') H 4> a (x - X) 



(3) 



and H is the approximate Hamiltonian when there is an 
applied field with vector potential A (x, t), then they are 
related by 



H(£ , ,£) = H (£',£)e 



iqA{t)-[X' -X)/hc 



with 



A(t) = (A(X',t)+A(X,t))/2. 



(4) 
(5) 



Here £ labels a localized basis function centered on a nu- 
cleus whose instantaneous position is X (£, t), and we 
adopt the convention of normally suppressing the indices 
I and £' , as well as the time t, by just writing X and 
X'. We will ignore any applied scalar potential Aq, any 
fis • B spin interactions, and the coupling of ion cores 
or nuclei to the applied fields, since these effects can be 
easily included when necessary. 

With the prescription of ([4]) , one does not need any new 
parameters in a calculation that employs either a semiem- 
pirical [II OJ] or first-principles-based [IS EE El El El 
Hamiltonian Hq whose elements are known as a function 
of (X — X'). On the other hand, this prescription is 
in one respect a rather crude approximation: It omits 
intra-atomic excitations, and would therefore give no ex- 
citation at all for isolated atoms. 

Here a more general version of the result of Ref. (2(| 
will be obtained, in a form which is almost equally conve- 
nient for large-scale applications, although it does require 
new parameters - namely dipole matrix elements 

Mo (£', £)=Q I d 3 x 0;, (x - X') (x - X) a (x - X) (6) 



and on-site (X' = X) matrix elements of the momentum 
operator 

Po {£', £) = ( d 3 x c/>* a , (x - X') (-ihV) cf> a (x - X) (7) 



where a labels an orbital centered on the nucleus whose 
instantaneous position is X. Recall that £ labels both 
nucleus and orbital, so at a given instant in time 

£ <-> X,a. (8) 

One key step is to expand ip in terms of London or- 
bitals, which we define to be any localized basis functions 
4> a that are related to field-independent basis functions 

<\>a by 



(f> a (x - X,t) 



a iqA{x,t)-(x—X) I he 



X) 



(9) 
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Notice that 4> a (x — X, t) = <f> a (x — X) when A = 0, so 
that after the application of a laser pulse, for example, 
the London orbitals return to being standard basis func- 
tions. The 4> a need not be a complete set, but should, of 
course, be a large enough set to model all physically rel- 
evant phenomena. The relatively weak time dependence 
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of the nuclear positions X is ignored for the moment, but 
will be included below. The original Hamiltonian of $Z$i 
can be rewritten as [13, H3 

jj _ e %q f A(x,t)-dx/hcjj^ e -iq f A(x,t)-dx/hc 



H = p 2 /2m + U , p=-ihV 



(11) 



since © and (flU)) yield the same result when operating 
on an arbitrary function, and are therefore the same op- 
erator. As will be seen immediately below, there are no 
problems in interpreting the integral of p0[) in the way 
that it is used here, since it is well-defined locally in eval- 
uating each matrix element. 

We now need the single approximation that A(x,t) 
varies slowly with respect to x over an atomic diameter 
or bond length, so that 



A(x,t) w A(t) 



(12) 



in the matrix elements which involve <f> a i (x — X') and 
4> a (x — X). (The wavelength is thus assumed to be large 
compared to 1 A.) When (TlT)|) and 



1p (X, t)=22<P V, t) <f>a (X - X, t) 
I 



(13) 



are substituted into |T]), and the resulting equation is 
subjected to / d 3 x <f)* a , (x — X' , t), we then obtain 



E 



S (£', £) ih d ^^; - = J2H {£', I) V> (I, t) (14) 



dt 



where 



S (£',£) 
So {£',£) 
H (£',£) 

fi (£',£) 



= S 



ft e iqA(t)-(x'-X)/hc 

Jd 3 x^* a , (x-X')<t> a {x 
H (£',£) -E(t) -n {£',£) 

Mo (l>J )e «lMt)iX'-X)/Hc 



X) 



and 



E(t) 



1 dA (t) 

c dt 



is the electric field. In matrix form, (14) is 



H-ij,{t). 



(15) 
(16) 

(17) 
(18) 

(19) 
(20) 



If there are N e electronic basis functions, then ip is an 
7V e -dimensional vector, whereas x, A, /x, etc. are 3- 
dimensional vectors. The dipole matrix elements can in 
principle be obtained in ab initio calculations like those 
used to obtain, e^g., the Hamiltonian matrix elements 
H (£', i) [!, i, 0, Sj. Alternatively, one might make the 
approximation of including only the terms with single- 
atom dipole matrix elements, (Xa' , Xa), and then 
take these from either atomic calculations or experiment. 



We now return to the time dependence of the nuclear 
positions X. With (Ti"3"|) rewritten as 

f/>(x,t) = ^2${l,t) <Pa(x-X) (21) 
t 

$(l,t) = t/> (£,t) e i9A(x,ty(x-X)/hc ( 22 ) 

we have [22[ 

'dip {£) 



dt 2-" 



dt 



■K{ X -x )+ Ht) d(t)Ax dx X) -x 



In order to treat the second term above, we assume 
(as indicated by the notation) that the basis functions 
depend only on (x — X), so that 



d<\> a (x - X) 
dX 



dcj> a [x - X) 
d{x-X) 
= -VJ> a {x-X). 



(23) 
(24) 



There is an additional correction involving X which 
arises from 



dip{£) 
dt 



^iqA(x)-(x — X) /he 



dip {£) 



dt 

dA (x) 
dt ' 



(x - X)-A(x) ■ X 



It follows that (fr?]) is modified to 

H (£',£) = H (e',£)e iqA i x '- x ) /hc 

-E(t)-n(£',£)-X-P(£',£) (25) 

where 

P (£',£) = p(l',£) + (q/c) AS (£',£) (26) 
p (£>,£) = p Q (£',£) e lqA i X '- X ) /hc (27) 

so another set of parameters is needed to treat the time 
dependence of the basis functions that arises from nuclear 
motion namely, the matrix elements of the momentum 
operator — ifiV. 

However, when X' =/= X, there is a more convenient 
way of writing p {£', £)■ 



po (£',£) 



ih J d 3 x(t)* a , (x - x') 



dX 



■(28) 



ih— S (£',£) if X'^X. (29) 



Furthermore, in the usual case of basis functions 
( "atomic orbitals" ) which are either even or odd under 
inversion through the nucleus, the fact that (x — X) and 
V = d/d (x — X) are odd under inversion (with X here 
taken to be fixed) implies that 



Mo (£,£)= po (£,£) = 0. 



(30) 
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Notice that (|25|) respects gauge invariance: If 

A (t) — » A' (t) — A(t) + AA (31) 
where AA is independent of t, then (fl4|) still holds with 



1>(£,t)->fl>' (£,t) 



_ iqAA-X/hc 



(32) 



This is the discrete version of 



A(x,t) -> A'(x,t) = A(x,t) + VA(x) (33) 
V- (sb, t) V' (as, *) = e iqA W Hc ip (x, t) . (34) 

If AA is a function of t, gauge invariance again holds, 
but with the scalar potential included. 

Equation (|25p is the central result of the present note. 
This effective Hamiltonian is not manifestly Hermitian, 
but it still conserves probability and preserves the Pauli 
principle, since a straightforward calculation using (|25j) 
in |T1|) gives 



m a (yl ■ s • v>„) /at = o 



(35) 



where n labels a time-dependent one-electron state. This 
result also follows from the original Schrodinger equation 
(Q]) and the expansion (fl"3)) . since 



d 3 xij*n>(x,t)i>n(x,t)=i>l(t)-S{t)-lP n (t) , (36) 



but it is reassuring that our approximation (|12p preserves 
orthonormality of the time-dependent states. 

For slowly moving nuclei the last term in (|25p is not 
important. (It may be worth mentioning in this context 
that the direct coupling of the nuclei to the field is not 
considered here, since it can be treated separately.) In an 
earlier paper [23j we argued that the nuclear motion can 
be approximately treated as a "nuclear velocity field" 
analogous to the radiation field, and in this spirit we 
obtained (as a crude approximation) a generalized Peierls 
substitution: 

H eff (t,£) = e MM(*')+™*T*' 

x H (e!,£)e-ili A W +m *]- x . (37) 

We also used this modified Hamiltonian in calculations 
for organic molecules responding to femtosecond-scale 
laser pulses of moderately strong intensity (<~ 10 12 
W/cm 2 ), and found that the X terms made very little 
difference in the final results. On the other hand, the 2- 
center momentum matrix elements can be obtained from 
(|29p . and the nonzero 1-center matrix elements from ei- 
ther atomic calculations or experiment, so it is certainly 
feasible to include the last term in (|25| . Notice that this 
term is different from the Pulay correction [l[ , which also 
results from the fact that the basis functions follow the 
nuclei, but occurs in the equation of motion for the nuclei 
rather than the time-dependent Schrodinger equation for 
the electrons. In the kind of approach considered here 



there is no Pulay correction, because the Hamiltonian 
matrix elements are supposed to have a position depen- 
dence that includes the movement of the basis functions. 

In this context, it is worth noting that the "Ehren- 
fest dynamics" [2J, [2f| of, e.g., time-dependent density- 
functional theory (TDDFT) and the density-functional- 
based calculations of Refs. [ll El El El El , 

can be sub- 
stantially improved in molecular calculations via a triv- 
ially different procedure that might be called "reduced 
Ehrenfcst dynamics" and which is similar in spirit to the 
surface hopping methods of Tully and others [II [27]]. Let 
us first recall some well-known results: The total wave- 
function for a system of nuclei, with coordinates X n , and 
electrons, with coordinates be represented by the 

Born-Oppenheimer expansion 

4- tot [X n , x e ,t) = J2 ^ i x n, t) {x e , X n ) . (38) 

i 

The basis functions ^ are the electronic eigenstates at 
fixed X n , with the electron- nuclei and nuclei- nuclei in- 
teractions included in the electronic Hamiltonian H e : 

H e [X^ (x e , X n ) = Ei [X n ) (x e ,X n ) . (39) 

Substitution into the Schrodinger equation 



ihdf tot /dt = HV* 



n = T n + H e , 



(40) 



where T n is the nuclear kinetic energy operator, gives an 
equation of the form [H, [2!| 



at 



n 2 



(41) 



2M n y ~ VJ 

(i\Vn\j 



{2F ij -V n + G ij ) (42) 
G« = (i|V»|j) (43) 



where M n is a representative nuclear mass and V n in- 
volves all the appropriately rescaled nuclear coordinates. 
If there are N n relevant nuclear coordinates, then V n 
and F are 7V„-dimensional vectors. Also, quantities in 
the last line are matrix elements defined in terms of 
and ^ j in the usual way. If the components $i are as- 
sembled into a vector <E>, (|4"Tj) can be written in a form 
which resembles a nonabelian gauge theory [3(3]: 



2M n 



(V„ 



E 



■ $ 



(44) 



where E is the diagonal matrix with elements Ei. Finally, 
it can be shown that [31j 



(*]Vnge|j) 

Ej — Ei 



Ei ? Ej 



(45) 



This last equation implies that each term in the Born- 
Oppenheimer expansion should evolve nearly indepen- 
dently if it is sufficiently distant in energy from all the 
other terms: If 



|^-^|>|(i|v n if e |i)|(V^) 



(46) 
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where P* = (2M n E l ) 1/2 , then 



(47) 



This is the time-dependent Born-Oppcnheimer or adia- 
batic approximation. 

On the other hand, whenever nuclear motion causes 
two Born-Oppcnhcimcr "potential energy surfaces" to 
approach each other, so that 



i^-£?ii<K<iv n H ii)i(w: 



(48) 



there is a nonadiabatic interaction [28|, |29j, |32|, |33j , and a 
Born-Oppenheimer simulation based on (1471) is invalid. 

The results of Refs. 0, El Ell EE Ejhave provided 
a clear demonstration of the following features of simu- 
lations based on Ehrenfest dynamics: 

(1) Electronic transitions are automatically observed 
at the points of closest approach where (|48|) holds, 
with energy released to molecular vibrations. These 
points are, of course, avoided crossings near the con- 
ical intersections in configuration space predicted by 
Teller [IE SI 13. 

(2) These transitions occur rapidly, over a time interval 
of ~ 1 femtosecond, during which the nuclei do not move 
appreciably. 

Ehrenfest simulations are based on the equation of mo- 
tion for the Heisenber g op erator X (t) representing any 
nuclear coordinate 



25J: 



Md 2 X/dt 2 = -dH/dX. 



(49) 



Here M is the corresponding nuclear mass and TL is the 
Hamiltonian of the system. In a standard Ehrenfest sim- 
ulation, the expectation value is taken over the full state 
of the system, including excitations (e.g. by a laser pulse) 
and de-excitations (e.g. by nuclear motion near conical 
intersections): 



M 



d 2 (X) 
dt 2 



mix 



dX 



on (<x>) 

d(X) 



(50) 



There are clearly two weaknesses with this approach: 
First, the equality on the left represents an average over 
all the terms in the expansion (|38p . with each term rep- 
resenting a different nuclear trajectory. Second, the ap- 
proximation on the right is totally invalid if these trajec- 
tories are very different. 

Suppose, however, that the standard procedure for an 
Ehrenfest simulation is replaced by a trivially different 
procedure in which the state of the system is collapsed 
to a single Born-Oppenheimer term immediately after 
an excitation or de-excitation event. Then (|47|) implies 
that it will essentially remain in this single adiabatically 
evolving state until the next such event. For this re- 
duced electronic state, the nuclei will ordinarily follow 
a single trajectory, except for quantum fluctuations of 

order ^ (^X — (X)^j ^ [35| . It is still possible for nu- 
clear wavepackets to diverge on a single potential energy 



surface, but one does not expect this to be a common oc- 
currence for processes in which the most relevant nuclei 
are reasonably heavy. 

For simplicity, first consider a very short laser pulse 
(e.g. ^1 — 5 femtoseconds in duration) applied to a 
molecule. The procedure for a reduced Ehrenfest simula- 
tion is as follows: Start with a single electronic eigenstate 
(e.g. the ground state) and initially perform an Ehrenfest 
simulation in the usual way. Immediately following the 
pulse, the molecule will be in a superposition of electronic 
eigenstates: 



(*) = £• 



(51) 



At this point one collapses \I/ e to a single eigenstate 
and continues the simulation, with (X) now interpreted 
as the expectation value for this single resulting time- 
dependent state, until another significant excitation or 
de-excitation is observed, after which there is again a 
further reduction to a single electronic eigenstate. 

There are potentially a substantial number of branches 
to be followed during this sort of simulation, correpond- 
ing to the various states in the superposition (|5ip after 
an excitation or de-excitation event. The goal, however, 
is to understand the most relevant processes, and there 
will ordinarily be physical motivations for selecting the 
most interesting branches. Similarly, there will be many 
branches emerging during an excitation process whose 
duration is long enough for the nuclei to move appre- 
ciably before it is completed (e.g., a femtosecond-scale 
laser pulse whose duration is still 3> 1 femtosecond) , and 
a choice among the branches again has to be based on 
physical considerations. 

For a molecule subjected to high-frequency or high- 
intensity radiation, the branches include ionized states. 
The one-electron matrix element between an orbital £ 
and an ionized state with momentum p is 



H, 



—A(x,t)-(e\p\p). 

mc 



(52) 



For a crude description of ionization, one might add a 
model orbital cf>o to the basis, with 



H t = — \A {X 7 t)\p Q 
mc 



H i0 = 



(53) 



where po ~ h/ao, ao is the Bohr radius, and c*i is an 
adjustable dimensionless parameter. This non-Hermitian 
Hamiltonian removes amplitude from the orbital £ at each 
time step and does not return it, so it crudely models 
excitation to a localized wavepacket with the electron 
ultimately escaping the system. An appreciable proba- 
bility for a given ionized state then provides motivation 
for following that branch in a reduced Ehrenfest simu- 
lation. Notice that an accurate treatment of ionization 
is not necessary if the only issue is whether an ionized 
state is important enough to warrant a simulation of the 
subsequent dynamics in that state. Also notice that the 
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energy Hqq of the extra orbital is irrelevant (so one can 
take H 00 = 0) and that a single extra orbital is sufficient 
regardless of the size of the system. 

After each wavefunction collapse, the use of ([50]) 
implies that the nuclei are treated classically. It is 
then approp riate to use the mixed classical-quantum ac- 
tion HE ii| S = fdtL, where 



j/ 4 n< 



h.c. 



-u, 



rep 



(54) 



where 7Y e is the electronic Hamiltonian, |\& e ) is the elec- 
tronic state, "h.c." means "Hermitian conjugate" , k la- 
bels a nucleus with spatial coordinates a, and U rep is 
the repulsive interaction between nuclei or ion cores. As 
shown in Ref. [36| (but with H now given by (|2"5")) ). ex- 
tremalization of this action leads to (j2"0|) and 



M ■ 



, rf 2 A 



1 (+ fdH ^dS d 



ax at 



-h.c. 



dX 



(55) 



if one makes the usual time-dependent effective-field ap- 
proximation, with exchange and correlation represented 
by an effective one-electron potential, and the electronic 
state represented by a single antisymmetrized product 
wavefunction ty e (t). Here X is any nuclear coordinate 
and M is the corresponding mass. 

The reduced Ehrenfest method described above com- 
bines the advantages of Born-Oppenhcimcr simulations, 
which are valid when (|4"T)l holds, and Ehrenfest simula- 
tions, which are suitable for treating the vibronic transi- 
tions when (|48|) holds, as the results of Refs. [HI, 

nam 

[HI, [l!| have clearly demonstrated. The use of reduced 



Ehrenfest simulations should solve various problems that 
are encountered in standard Ehrenfest simulations - for 
example, the apparent failure of TDDFT to correctly 
describe the isomerization of retinal Q. One problem 
with TDDFT is that the energies of excited states are 
not accurately described, but a potentially more severe 
problem in the case of molecules is that TDDFT is a 
special case of standard Ehrenfest dynamics, and as a re- 
sult fails to yield a complete return to the ground state 
following de-excitation near a conical intersection. In a 
reduced Ehrenfest simulation, on the other hand, one 
correctly follows the nuclear dynamics for that fraction 
of the population of molecules which does return to the 
ground state, and which therefore should isomerize more 
readily. Reduced Ehrenfest simulations are practical for 
large molecules, and are still consistent with the true 
meaning of quantum amplitudes, which yield probabil- 
ities for the various outcomes that are observed at the 
classical level. 

Finally, it may be worth noting that the above treat- 
ment can be straightforwardly generalized to other par- 
ticles, relativistic systems, and nonabelian gauge fields, 
with ip in (TT]) interpreted as a multicomponent field and 
the Hamiltonian of ([2]) appropriately changed. It can also 
be used with many-body effects included through self- 
energy terms, in the Kadanoff-Baym/Keldysh e quat ions 
for time-dependent and nonequilibrium problems |37l.l38l|. 
The chief limitation is the use of localized basis functions 
and the approximation (|12p . 
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